module operators_mod

	use const_mod
	use mesh_mod
	use state_mod
	use tracer_mod
	use parallel_mod
	use vert_coord_mod

implicit none

	private
	public calc_ph
	public calc_m
	public calc_mf
	public calc_hori_uv
	public calc_hori_div
	public calc_hori_cfl
	public calc_vert_cfl
	public calc_vert_we
	public calc_we_split_wgt
	public calc_we_split_wgt2

contains	

	subroutine calc_ph(state)

		type(state_type), intent(inout) :: state

		integer i, j, k

		associate (phs    => state%phs   , & !in
							 ph_lev => state%ph_lev, & !out
							 ph     => state%ph    )   !out

			do k = mesh%half_lev_ibeg, mesh%half_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						ph_lev(i,j,k) = hybrid_coord_calc_ph_lev(k, phs(i,j))
					end do
				end do
			end do

			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						ph(i,j,k) = hybrid_coord_calc_ph(k, phs(i,j))
					end do
				end do
			end do
	
		end associate

	end subroutine calc_ph

	subroutine calc_m(state)

		type(state_type), intent(inout) :: state

		integer i,j,k

		associate (ph_lev => state%ph_lev, & !in
			              m => state%m     , & !out
			          m_lon => state%m_lon , & !out
			          m_lat => state%m_lat , & !out
			          m_lev => state%m_lev)
			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						m(i,j,k) = ph_lev(i,j,k+1) - ph_lev(i,j,k)
					end do
				end do
			end do
			call fill_zonal_halo_cell(m, all_halo=.true.)

			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%half_lon_ibeg, mesh%half_lon_iend
						m_lon(i,j,k) = 0.5 * (m(i,j,k) + m(i+1,j,k))
					end do
				end do
			end do
			call fill_zonal_halo_lon(m_lon, all_halo=.true.)

			do k = mesh%full_lev_ibeg, mesh%full_lev_iend
				do j = mesh%half_lat_ibeg, mesh%half_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						m_lat(i,j,k) = 0.5 * (m(i,j,k) + m(i,j+1,k))
					end do
				end do
			end do
			call fill_zonal_halo_lat(m_lat, all_halo=.true.)

			do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend -1
				do j = mesh%full_lat_ibeg, mesh%full_lat_iend
					do i = mesh%full_lon_ibeg, mesh%full_lon_iend
						m_lev(i,j,k) = 0.5 * (m(i,j,k) + m(i,j,k-1))
					end do
				end do
			end do
			m_lev(:,:,mesh%half_lev_ibeg) = 0.5 * m(:,:,1)
			m_lev(:,:,mesh%half_lev_iend) = 0.5 * m(:,:,mesh%full_lev_iend)

		end associate

	end subroutine calc_m

	subroutine calc_mf(state)

		type(state_type), intent(inout) :: state

		integer i, j, k

		associate (u => state%u,     & !in
			         v => state%v,     & !in
			     m_lon => state%m_lon, & !in
			     m_lat => state%m_lat, & !in
			    mf_lon => state%mf_lon,& !out
			    mf_lat => state%mf_lat)  !out
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			do j = mesh%full_lat_ibeg, mesh%full_lat_iend
				do i = mesh%half_lon_ibeg, mesh%half_lon_iend
					mf_lon(i,j,k) = m_lon(i,j,k) * u(i,j,k)
				end do
			end do
		end do
		call fill_zonal_halo_lat(mf_lon, all_halo=.true.)
		
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
			do j = mesh%half_lat_ibeg, mesh%half_lat_iend
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
					mf_lat(i,j,k) = m_lat(i,j,k) * v(i,j,k)
				end do
			end do
		end do
		end associate

	end subroutine calc_mf

  subroutine adv_accum_mf_hori(old_state, new_state, old_tracer, dt)

  	type(state_type), intent(in) :: old_state
  	type(state_type), intent(in) :: new_state
  	type(tracer_type), intent(inout) :: old_tracer
  	real(8), intent(in) :: dt
  	integer i, j, k

  	associate (mfx => old_tracer%mfx, mfy => old_tracer%mfy)
  		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
  			do j = mesh%full_lat_ibeg, mesh%full_lat_iend
  				do i = mesh%half_lon_ibeg, mesh%half_lon_iend
  					mfx(i,j,k) = old_state%mf_lon(i,j,k) + new_state%u(i,j,k) * new_state%m_lon(i,j,k)
  				end do
  			end do
  		end do

  		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
  			do j = mesh%half_lat_ibeg, mesh%half_lat_iend
  				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  					mfy(i,j,k) = old_state%mf_lat(i,j,k) + new_state%v(i,j,k) * new_state%m_lat(i,j,k)
  				end do
  			end do
  		end do

  	end associate

  end subroutine adv_accum_mf_hori

	subroutine calc_hori_div(old_tracer, dt)

		type(tracer_type), intent(inout) :: old_tracer
		real(8), intent(in) :: dt
		real(r8) pole
		integer i, j, k

		associate (u    => old_tracer%u   , &
			 				 v    => old_tracer%v   , &
							 divx => old_tracer%divx, &
							 divy => old_tracer%divy)

	  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
	  	do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
	  		do i = mesh%full_lon_ibeg, mesh%full_lon_iend
	  			divx(i,j,k) = (u(i,j,k) - u(i-1,j,k)) / mesh%de_lon(j)
	  			if (abs(divx(i,j,k) * dt) > 1) then
	  				print*, abs(divx(i,j,k) * dt), 'The zonal deformational Courant number exceeds 1'
	  			end if
	  			divy(i,j,k) = (v(i,j  ,k) * mesh%le_lat(j  ) - &
	  									   v(i,j-1,k) * mesh%le_lat(j-1)) / mesh%area_cell(j)
	  			if (abs(divy(i,j,k) * dt) > 1) then
	  				print*, abs(divy(i,j,k) * dt), 'The meridional deformational Courant number exceeds 1'
	  			end if
	  		end do
	  	end do
    end do
  	!south pole
  	do k = mesh%full_lev_ibeg, mesh%full_lev_iend
  	  j = mesh%full_lat_ibeg
  	  pole = sum(v(mesh%full_lon_ibeg:mesh%full_lon_iend,j,k)) * mesh%le_lat(j) / mesh%num_full_lon / mesh%area_cell(j)
  	  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  	  	divy(i,j,k) = pole
  	  end do
  	end do

  	! north pole
  	do k = mesh%full_lev_ibeg, mesh%full_lev_iend
	  	j = mesh%full_lat_iend
	  	pole = sum(v(mesh%full_lon_ibeg:mesh%full_lon_iend,j-1,k)) * mesh%le_lat(j-1) / mesh%num_full_lon / mesh%area_cell(j)
	  	do i = mesh%full_lon_ibeg, mesh%full_lon_iend
	  		divy(i,j,k) = -pole
	  	end do
	  end do
  	end associate 

	end subroutine calc_hori_div
	
	subroutine calc_hori_uv(old_tracer, old_state, new_state, dt)
		type(tracer_type), intent(inout) :: old_tracer
		type(state_type), intent(in) :: old_state
		type(state_type), intent(in) :: new_state
		real(8), intent(in) :: dt 
		integer i, j, k

		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
  	  do j = mesh%full_lat_ibeg, mesh%full_lat_iend
  	  	do i = mesh%half_lon_ibeg, mesh%half_lon_iend
  	  		old_tracer%u(i,j,k) = (old_state%u(i,j,k) + new_state%u(i,j,k)) * 0.5_r8
  	  	end do
  	  end do
  	end do
    call fill_zonal_halo_lon(old_tracer%u, all_halo=.true.)
		
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
  	  do j = mesh%half_lat_ibeg, mesh%half_lat_iend
  	  	do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  	  		old_tracer%v(i,j,k) = (old_state%v(i,j,k) + new_state%v(i,j,k)) * 0.5_r8
  	  	end do
  	  end do
  	end do
  	call fill_zonal_halo_lat(old_tracer%v, all_halo=.true.)

	end subroutine calc_hori_uv

	subroutine calc_hori_cfl(old_tracer, dt)

		type(tracer_type), intent(inout) :: old_tracer
		real(8), intent(in) :: dt 
		integer i, j, k

  	associate (u    => old_tracer%u   , &
  						 v    => old_tracer%v   , &
  						 cflx => old_tracer%cflx, &
  	  				 cfly => old_tracer%cfly)
		do k = mesh%full_lev_ibeg, mesh%full_lev_iend
	    do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
	    	do i = mesh%half_lon_ibeg, mesh%half_lon_iend
	    		cflx(i,j,k) = u(i,j,k) * dt / mesh%de_lon(j)
	    		if (abs(cflx(i,j,k)) > mesh%lon_halo_width) then
	    	    print*, 'maxval of cflx', maxval(abs(cflx)), 'larger than lon_halo_width', mesh%lon_halo_width
	    	    stop
	        end if
	    	end do
	    end do
	  end do

	  do k = mesh%full_lev_ibeg, mesh%full_lev_iend
	    do j = mesh%half_lat_ibeg, mesh%half_lat_iend
	    	do i = mesh%full_lon_ibeg, mesh%full_lon_iend
	    		cfly(i,j,k) = v(i,j,k) * dt / mesh%de_lat(j)
	    	end do
	    end do
	  end do
	  if (maxval(abs(cfly)) > 1) then
	  	print*, 'cfly', maxval(abs(cfly)), 'exceeds 1!'
	  	stop
	  end if
	  end associate

	end subroutine calc_hori_cfl

  subroutine calc_vert_cfl(we, m_lev, cflz, dt)

  	real(r8), intent(in) :: we(mesh%full_lon_lb:mesh%full_lon_ub, &
													     mesh%full_lat_lb:mesh%full_lat_ub, &
													     mesh%half_lev_lb:mesh%half_lev_ub)
  	real(r8), intent(in) :: m_lev(mesh%full_lon_lb:mesh%full_lon_ub, &
											   		      mesh%full_lat_lb:mesh%full_lat_ub, &
															    mesh%half_lev_lb:mesh%half_lev_ub)
  	real(r8), intent(out) :: cflz(mesh%full_lon_lb:mesh%full_lon_ub, &
											   		      mesh%full_lat_lb:mesh%full_lat_ub, &
															    mesh%half_lev_lb:mesh%half_lev_ub)
  	real(8), intent(in) :: dt
  	integer i, j, k

	  do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
	  	do j = mesh%full_lat_ibeg, mesh%full_lat_iend
	  		do i = mesh%full_lon_ibeg, mesh%full_lon_iend
	  			if (we(i,j,k) > 0) then
	  				cflz(i,j,k) = we(i,j,k) * dt / m_lev(i,j,k) / (mesh%full_dlev(k-1) / mesh%half_dlev(k))
	  			else if (we(i,j,k) < 0) then
	  				cflz(i,j,k) = we(i,j,k) * dt / m_lev(i,j,k) / (mesh%full_dlev(k) / mesh%half_dlev(k))
	  			else
	  				cflz(i,j,k) = 0
	  			end if
	  		end do
	  	end do
	  	! if (maxval(abs(cflz(:,:,k))) > 1) then
	  	! 	print*, 'cflz=', maxval(abs(cflz(:,:,k))), ' exceeds 1!'
	  	! end if
	  end do

  end subroutine calc_vert_cfl

 	subroutine calc_vert_we(old_tracer, old_state, new_state)

 		type(tracer_type), intent(inout) :: old_tracer
 		type(state_type), intent(in) :: old_state
 		type(state_type), intent(in) :: new_state
 		integer i, j, k
 		
 		do k = mesh%half_lev_ibeg, mesh%half_lev_iend
 			do j = mesh%full_lat_ibeg, mesh%full_lat_iend
 				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
 					old_tracer%we(i,j,k) = 0.5_r8 * (old_state%we(i,j,k) + new_state%we(i,j,k))
 				end do
 			end do
 		end do

 	end subroutine calc_vert_we

  subroutine calc_we_split_wgt(old_tracer, old_state, dt)

  	type(tracer_type), intent(inout) :: old_tracer
  	type(state_type), intent(in) :: old_state
  	real(8), intent(in) :: dt
  	real(r8), parameter :: alpha_max = 1.1, alpha_min = 0.8, kesi = 0.9
  	real(r8) alpha_h, alpha_v, alpha_star_max, alpha_star_min
  	real(r8) pole, beta
  	integer i, j, k

  	associate (u      => old_tracer%u     , & ! in 
  						 v      => old_tracer%v     , & ! in
  						 we     => old_tracer%we    , & ! in
  						 m_lev  => old_state%m_lev  , & ! in
  						 we_exp => old_tracer%we_exp, & ! out
  						 we_imp => old_tracer%we_imp)   ! out

  	do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
  		do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				if (we(i,j,k) >= 0.0) then
          	alpha_h = dt * ((max(u(i,j,k-1), 0.) - min(u(i-1,j,k-1), 0.)) * mesh%le_lon(j) + &
  											    (max(v(i,j,k-1), 0.) * mesh%le_lat(j) - min(v(i,j-1,k-1), 0.) * mesh%le_lat(j-1))) / mesh%area_cell(j)
  				else
						alpha_h = dt * ((max(u(i,j,k), 0.) - min(u(i-1,j,k), 0.)) * mesh%le_lon(j) + &
  											    (max(v(i,j,k), 0.) * mesh%le_lat(j) - min(v(i,j-1,k), 0.) * mesh%le_lat(j-1))) / mesh%area_cell(j)
  				end if
  				alpha_star_max = alpha_max - kesi * alpha_h
  				alpha_star_min = alpha_min * alpha_star_max / alpha_max
  				alpha_v = dt * abs(we(i,j,k) / m_lev(i,j,k))
  				if (alpha_v <= alpha_star_min) then
  					beta = 1.0_r8
  				else if (alpha_v > alpha_star_min .and. alpha_v <= 2.0 * alpha_star_max - alpha_star_min) then
  					beta = 1.0_r8 / (1.0 + (alpha_v - alpha_star_min)**2 / (4.0 * alpha_star_max * (alpha_star_max - alpha_star_min)))
  				else
  					beta = alpha_star_max / alpha_v
  				end if
  				if (beta < 0 .or. beta > 1) then
  					print*, "vertical split weight out range from 0 to 1!"
  					stop
  				end if
  				we_exp(i,j,k) = beta * we(i,j,k)
  				we_imp(i,j,k) = (1.0 - beta) * we(i,j,k)
  			end do
  		end do
  	end do

  	!south pole
  	do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
			j = mesh%full_lat_ibeg
			pole = 0
			if (we(mesh%full_lon_ibeg,j,k) >= 0) then
			  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
			  	pole = pole + max(v(i,j,k-1), 0.)
			  end do
			else
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
			  	pole = pole + max(v(i,j,k), 0.)
			  end do
			end if
			alpha_h = dt * pole * mesh%le_lat(j) / mesh%num_full_lon / mesh%area_cell(j)
			alpha_star_max = alpha_max - kesi * alpha_h
			alpha_star_min = alpha_min * alpha_star_max / alpha_max
			alpha_v = dt * abs(we(mesh%full_lon_ibeg,j,k) / m_lev(mesh%full_lon_ibeg,j,k))
			if (alpha_v <= alpha_star_min) then
				beta = 1.0
			else if (alpha_v > alpha_star_min .and. alpha_v <= 2.0 * alpha_star_max - alpha_star_min) then
				beta = 1.0 / (1.0 + (alpha_v - alpha_star_min)**2 / (4.0 * alpha_star_max * (alpha_star_max - alpha_star_min)))
			else
				beta = alpha_star_max / alpha_v
			end if
	  	if (beta < 0 .or. beta > 1) then
				print*, "vertical split weight out range from 0 to 1!"
				stop
			end if
			we_exp(:,j,k) = beta * we(mesh%full_lon_ibeg,j,k)
			we_imp(:,j,k) = (1.0 - beta) * we(mesh%full_lon_ibeg,j,k)
		end do

  	!north pole
  	do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
  		j = mesh%full_lat_iend
  		pole = 0
  		if (we(mesh%full_lon_ibeg,j,k) >= 0) then
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				pole = pole - min(v(i,j,k-1), 0.)
  			end do
  		else
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				pole = pole - min(v(i,j,k), 0.)
  			end do
  		end if
  		alpha_h = dt * pole * mesh%le_lat(j-1) / mesh%num_full_lon / mesh%area_cell(j)
  		alpha_star_max = alpha_max - kesi * alpha_h
  		alpha_star_min = alpha_min * alpha_star_max / alpha_max
  		alpha_v = dt * abs(we(mesh%full_lon_ibeg,j,k) / m_lev(mesh%full_lon_ibeg,j,k))
			if (alpha_v <= alpha_star_min) then
				beta = 1.0
			else if (alpha_v > alpha_star_min .and. alpha_v <= 2.0 * alpha_star_max - alpha_star_min) then
				beta = 1.0 / (1.0 + (alpha_v - alpha_star_min)**2 / (4.0 * alpha_star_max * (alpha_star_max - alpha_star_min)))
			else
				beta = alpha_star_max / alpha_v
			end if
	  	if (beta < 0 .or. beta > 1) then
				print*, "vertical split weight out range from 0 to 1!"
				stop
			end if
			we_exp(:,j,k) = beta * we(mesh%full_lon_ibeg,j,k)
			we_imp(:,j,k) = (1.0 - beta) * we(mesh%full_lon_ibeg,j,k)
		end do
		we_exp(:,:,mesh%half_lev_ibeg) = 0.0; we_exp(:,:,mesh%half_lev_iend) = 0.0
		we_imp(:,:,mesh%half_lev_ibeg) = 0.0; we_imp(:,:,mesh%half_lev_iend) = 0.0
		end associate
		
  end subroutine calc_we_split_wgt

  subroutine calc_we_split_wgt2(old_tracer, old_state, dt) 
    ! maybe has bug, fix me
  	type(tracer_type), intent(inout) :: old_tracer
  	type(state_type), intent(in) :: old_state
  	real(8), intent(in) :: dt
  	real(r8), parameter :: alpha_max = 1.1, alpha_min = 0.8, kesi = 0.9
  	real(r8) alpha_h, alpha_v, alpha_star_max, alpha_star_min
  	real(r8) pole, beta, dh
  	integer i, j, k

  	associate (u      => old_tracer%u     , & ! in 
  						 v      => old_tracer%v     , & ! in
  						 we     => old_tracer%we    , & ! in
  						 m_lev  => old_state%m_lev  , & ! in
  						 we_exp => old_tracer%we_exp, & ! out
  						 we_imp => old_tracer%we_imp)   ! out

  	do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
  		do j = mesh%full_lat_ibeg + 1, mesh%full_lat_iend - 1
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				if (we(i,j,k) >= 0.0) then
          	alpha_h = (max(u(i,j,k-1), 0.) - min(u(i-1,j,k-1), 0.)) * mesh%le_lon(j) * mesh%full_dlev(k-1) + &
  										(max(v(i,j,k-1), 0.) * mesh%le_lat(j) - min(v(i,j-1,k-1), 0.) * mesh%le_lat(j-1)) * mesh%full_dlev(k-1)
  					dh  = mesh%full_dlev(k-1)
  				else
						alpha_h = (max(u(i,j,k), 0.) - min(u(i-1,j,k), 0.)) * mesh%le_lon(j) * mesh%full_dlev(k) + &
  										(max(v(i,j,k), 0.) * mesh%le_lat(j) - min(v(i,j-1,k), 0.) * mesh%le_lat(j-1)) * mesh%full_dlev(k)
  					dh = mesh%full_dlev(k)
  				end if
  				alpha_star_max = alpha_max * dh - alpha_h * kesi * dt / mesh%area_cell(j)
  				if (alpha_star_max > 0) then
  					alpha_star_min = alpha_min * alpha_star_max / alpha_max
  				  alpha_v = abs(we(i,j,k)) * kesi * dt
  				  if (alpha_v < alpha_star_min) then
  					  beta = alpha_star_max**2
  				  else if (alpha_v < (2.0 - alpha_min / alpha_max) * alpha_star_max) then
  					  beta = alpha_star_max**2 + (1.0_r8 / (4.0 - 4.0 * alpha_min / alpha_max)) * (alpha_v - alpha_star_min)**2
  					else
  						beta = alpha_star_max * alpha_v
  					end if
  					we_exp(i,j,k) = alpha_star_max**2 * we(i,j,k) / beta
  					we_imp(i,j,k) = we(i,j,k) - we_exp(i,j,k)
  				else
  					we_exp(i,j,k) = 0.0
  					we_imp(i,j,k) = we(i,j,k)
  				end if
  			end do
  		end do
  	end do

  	!south pole
  	do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
			j = mesh%full_lat_ibeg
			pole = 0
			if (we(mesh%full_lon_ibeg,j,k) >= 0) then
			  do i = mesh%full_lon_ibeg, mesh%full_lon_iend
			  	pole = pole + max(v(i,j,k-1), 0.) * mesh%le_lat(j) * mesh%full_dlev(k-1)
			  end do
			  dh = mesh%full_dlev(k-1)
			else
				do i = mesh%full_lon_ibeg, mesh%full_lon_iend
			  	pole = pole + max(v(i,j,k), 0.) * mesh%le_lat(j) * mesh%full_dlev(k)
			  end do
			  dh = mesh%full_dlev(k)
			end if
			alpha_h = pole 
			alpha_star_max = alpha_max * dh - alpha_h * kesi * dt / mesh%num_full_lon / mesh%area_cell(j)
			if (alpha_star_max > 0) then
				alpha_star_min = alpha_min * alpha_star_max / alpha_max
			  alpha_v = abs(we(mesh%full_lon_ibeg,j,k)) * kesi * dt
			  if (alpha_v < alpha_star_min) then
				  beta = alpha_star_max**2
			  else if (alpha_v < (2.0 - alpha_min / alpha_max) * alpha_star_max) then
				  beta = alpha_star_max**2 + (1.0_r8 / (4.0 - 4.0 * alpha_min / alpha_max)) * (alpha_v - alpha_star_min)**2
				else
					beta = alpha_star_max * alpha_v
				end if
			  we_exp(:,j,k) = alpha_star_max**2 * we(mesh%full_lon_ibeg,j,k) / beta
  			we_imp(:,j,k) = we(mesh%full_lat_ibeg,j,k) - we_exp(mesh%full_lon_ibeg,j,k)
			else
				we_exp(:,j,k) = 0.0
				we_imp(:,j,k) = we(mesh%full_lon_ibeg,j,k)
			end if
		end do

  	!north pole
  	do k = mesh%half_lev_ibeg + 1, mesh%half_lev_iend - 1
  		j = mesh%full_lat_iend
  		pole = 0
  		if (we(mesh%full_lon_ibeg,j,k) >= 0) then
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				pole = pole - min(v(i,j,k-1), 0.) * mesh%le_lat(j-1) * mesh%full_dlev(k-1)
  			end do
  			dh = mesh%full_dlev(k-1)
  		else
  			do i = mesh%full_lon_ibeg, mesh%full_lon_iend
  				pole = pole - min(v(i,j,k), 0.) * mesh%le_lat(j-1) * mesh%full_dlev(k)
  			end do
  			dh = mesh%full_dlev(k)
  		end if
  		alpha_h = pole
   		alpha_star_max = alpha_max * dh - alpha_h * kesi * dt / mesh%num_full_lon / mesh%area_cell(j)
  		if (alpha_star_max > 0) then
  		  alpha_star_min = alpha_min * alpha_star_max / alpha_max
			  alpha_v = abs(we(mesh%full_lon_ibeg,j,k)) * kesi * dt
				if (alpha_v < alpha_star_min) then
				  beta = alpha_star_max**2
			  else if (alpha_v < (2.0 - alpha_min / alpha_max) * alpha_star_max) then
				  beta = alpha_star_max**2 + (1.0_r8 / (4.0 - 4.0 * alpha_min / alpha_max)) * (alpha_v - alpha_star_min)**2
				else
					beta = alpha_star_max * alpha_v
				end if
			  we_exp(:,j,k) = alpha_star_max**2 * we(mesh%full_lon_ibeg,j,k) / beta
  			we_imp(:,j,k) = we(mesh%full_lon_ibeg,j,k) - we_exp(mesh%full_lon_ibeg,j,k)
  		else
  		  we_exp(:,j,k) = 0.0
  			we_imp(:,j,k) = we(mesh%full_lon_ibeg,j,k)
  		end if
		end do
		we_exp(:,:,mesh%half_lev_ibeg) = 0.0; we_exp(:,:,mesh%half_lev_iend) = 0.0
		we_imp(:,:,mesh%half_lev_ibeg) = 0.0; we_imp(:,:,mesh%half_lev_iend) = 0.0
		end associate
		
  end subroutine calc_we_split_wgt2

end module operators_mod